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Abstract 

The concept of “numerical speed of 
so u n (L* is proposed in the construction of nu- 
merical flux. It is shown that this variable 
is responsible for the accurate resolution of 
discontinuities, such as contacts and shocks. 
Moreover, this concept can be readily ex- 
tended to deal with low speed and multiphase 
flows. As a result, the numerical dissipa- 
tion for low speed flows is scaled with the lo- 
cal fluid speed, rather than the sound speed. 
Hence, the accuracy is enhanced, the correct 
solution recovered, and the convergence rate 
improved. We also emphasize the role of mass 
flux and analyze the behavior of this flux. 
Study of mass flux is important because the 
numerical diffusivity introduced in it can be 
identified. In addition, it is the term common 
to all conservation equations. 

We show calculated results for a wide va- 
riety of flows to validate the effectiveness of 
using the numerical speed of sound concept 
in constructing the numerical flux. We espe- 
cially aim at achieving these two goals: (1) 
improving accuracy and (2) gaining conver- 
gence rates for all speed ranges. We find that, 
while the performance at high speed range is 
maintained, the flux now has the capability 


of performing well even with the low speed 
flows. Thanks to the new numerical speed of 
sound, the convergence is even enhanced for 
the flows outside of the low speed range. To 
realize the usefulness of the proposed method 
in engineering problems, we have also per- 
formed calculations for complex 31) turbulent 
flows and the results are in excellent agree- 
ment with data. 

1. Introduction 

Numerical representation of inviscid 
fluxes, namely the numerical flux function, 
has been a sub ject of intensive effort by many 
researchers during the last three decades. 
Despite the enormous progress that has been 
achieved, further analysis and deeper un- 
derstanding into these numerical procedures 
continue to draw interest and the findings 
are being reported. In this paper we will in- 
troduce the concept of the numerical speed 
of sound. This concept t urns out to be very 
useful in designing a numerical flux in order 
to meet certain criteria. Moreover, as will be 
seen later, the numerical speed of sound lends 
itself nicely to the formulation of numerical 
schemes for all speeds. In other words, it 
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plays an important role in bridging the gap 
between discretization schemes suitable for 
incompressible flows and those suitable for 
compressible flows. 


to be an important foundation for the devel- 
opment of recent schemes such as the AUSM- 
family [2-4], The mass flux of the Van Leer- 
type scheme is. 


For illustrative purposes, we shall begin 
by considering the ID flux for ideal gas. The 
in viscid numerical flux function is written as 
a sum of convective and pressure fluxes: 


f = pll 



= m 



It is noted that a, common mass flux m 
appears in all equations. This is also true 
for multidimensions. Since the mass flux is 
common for all equations, its effects will thus 
perpetuate in all variables. Hence, we suggest 
that it is desirable to observe this fact at the 
discrete level when devising a new scheme. 
However, this fact is not entirely enforced in 
several modern numerical schemes. 

The AUSMTamily schemes writes the in- 
terface flux ft/ 2 * in a form mimicking the con- 
tinuum flux, Eq. (1). as 


1 ' 

\ 

° ) 

a 

j + 

Pi/2 

Uy 

j/j+i 

0 / 


where the cell interface straddles the cells j 
and j + 1 and the subscript j/j + 1 denotes 
that the vector (1. t/. H) is evaluated with ei- 
ther j or j + 1 values according to whether 
ni i /2 is positive or negative. The detailed de- 
scription of m i/9 and pi/ 2 can be found in 
[2-4], 


2. Numerical Speed of Sound 


To understand the idea of numerical 
speed of sound, we firsl consider the cele- 
brated Van Leer scheme [1], which turns out 


1,1 i/2— Pj a jM^ 4 d) { Mj )+pj+\a J+ iM( 4 i i )) (Mj+ 1 ). 

( 3 ) 

where ^ are the split Mach number func- 
tions defined by the 4th-degree (or second de- 
gree in [1]) polynomials in the subsonic range: 


■K.JM) = 


Mf„. if \M\ > 1. 


where 

M± ) = l -(M±\M\). (5) 

and 

M% = ± ] -(M±l) 2 .^ (6) 

Other forms are possible, but are not es- 
sential in the discussion in this paper. The 
scheme. Eq. (3), is known to produce a dis- 
sipative result at a contact discontinuity. As 
M — M, — 4/,+i = 0. the mass flux does not 
vanish. 


"*1/2= ~ Pj+i (l J+i I/ 1 - (7) 

It is clear in the Van Leer scheme that two 
variables, namely (p, a), need to be dealt with 
if a mass flux is t o vanish at M = 0. It is due 
to this observation that a common speed of 
sound (i\j 2 is introduced in the Al SM + and 
AUSMDV schemes. Let s modify Eq. (3) by 
using « 1/2 . 

n?i/2= (i\ fi)( A/j )Tpj+i (.j i ))• 

. (S) 

where the exact definition of a 1/2 is not im- 
portant for the moment, but is certainly the 
central topic of this paper. The limiting form 
of nil / 2 at. M — 0 becomes. 

m l/2= a 1 / 2 ( Pj — Pj+ 1 )/ 1 0- (9) 


§The coefficient in (19b) of JCP 129 , 364-382 (1996), should be 1/4, not 1/2. 
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Clearly t he common speed of sound a \/- 2 alone 
is not enough to make the mass flux vanish 
at M = 0, and we need something extra. 
A novel approach is proposed in the AUS- 
MDV scheme in which the split Mach number 
functions M ± ^ are added as an anti-diffusive 
mechanism to fully cancel the excessive dis- 
sipation. The interested reader is referred to 
[4] for details. 

In a different approach, AUSM+ defines 
the mass flux as 


m 


1/2 = 


h/2 


d / 1 /2 ( Pj + Pj+ i ) “ | ^A/2 1 ( Pj+ 1 —Pj )] - 

( 10 ) 


where 


M \/2 — ). ( 11 ) 

Since XI \/2 — 0 when Mj = Mj+i — 0, we 
automatically get = 0. 


'2 

i o !«•••« 

8 - 
6 




0.8 
0.4 
G.O ' 
- 0.4 
-C.8 


Van ^eer 
HLlL 


4 8 

X 


dissipative schemes. Van Leer and HLLE [5], 
will destroy the discontinuity at the first time 
step: Van Leer flux is seen to be more dissi- 
pative than the HLLE flux. 

While a stationary contact can be exact ly 
captured by the original AUSM scheme [2], 
failure is encountered in the case of mov- 
ing contact discontinuity, as shown in Fig. 
2. This is also the reason for introduc- 
ing a common speed of sound in the AUSM 
scheme to achieve this and other improve- 
ments, thus leading to the improved scheme 
(ailed A1SM + [4]. Since the common speed 
of sound is int roduced as a means to yield an 
accurate solution, which is only meaningful 
in the numerical framework, hence it is here- 
after termed "numerical" speed of sound. 
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Fig. 2 Slowly moving contact; comparison of 
AUSM and AUSM+ solutions. 


Fig. 1 Stationary contact discontinuity by the 
Van Leer and HLLE schemes 

We remark that both AUSMDV and 
A I SM + yield the desired result for a station- 
ary contact discontinuity of any strength in 
(Pj'Pj+ 1) with an arbitrary ai/- 2 . Figure 1 dis- 
plays the solution at / > 0 for a stationary 
contact discontinuity. It is seen that the two 


As stated earlier, the contact (station- 
ary and moving) discontinuities can be ac- 
curately resolved with an arbitrary numerical 
speed of sound a j/ 2 . This leaves us one degree 
of freedom to determine a x / 2 so that another 
interesting property can he met, for example, 
exact, capturing of a stationary shock discon- 
tinuity. This is accomplished in the AUSM- 
family schemes by setting (detailed deriva- 
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lion can be found in [3]), 

(i \ ji — min(«i. «■/?), (7 = f/* 2 /niax(</*, |«|), 

( 12 ) 

where «* is t lie critical speed of sound. For 
ideal gas, we have 


*2 

a 


2(7-1)// 
1 + 1 


(13) 
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made to satisfy the same property by using 
the formula in Eq. (12) as well. Thirdly, in 
addition to the familiar Roe averaged speed 
of sound, we find that 


a 1/2 = max(fl L ,flfl.). (1-1) 

instead of "min" used in Eq. (12), can also 
yield an exact shock profile. The result in 
Fig. 3 for the Roe scheme was obtained with 
this formula. To our knowledge, these choices 
of a 1/2 to achieve the exact property with the 
HLLE and Roe schemes are previously un- 
known in the literature. 

Formulation for All Speeds 
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Fig. 3 Stationary shock discontinuity: a numer- 
ical speed of sound is inserted in each of the 
schemes included. 

It is well known that a single shock can 
be captured by the Roe scheme with the so- 
called Roe averaged quantities. One of these 
(iuaut.il ies is the averaged speed of sound, 
which also is "numerical*’ in nature. The 
idea of regarding the quantity a^/ 2 as a free 
variable turns out to be rather powerful in 
the sense that the idea can be applied to 
other schemes equally well. Firstly, the two 
intermediate shock points by the Van Leer 
scheme, including its original form [1] and 
H artel's modification [6] for the energy flux, 
are now dramatically removed by simply in- 
corporating this special numerical speed of 
sound. Secondly, the HLLE schemes can he 


It is widely known that the standard form 
of compressible equations, discretized with ei- 
ther centered or upwind schemes, yields two 
major effect s on the solution as the flow speed 
approaches zero: (1) a drastic slowdown or 
level-off of convergence rate, (2) an inaccu- 
rate or even incorrect solution. An effective 
way of dealing with the first point is by in- 
serting a premultiplying matrix to the time- 
derivative term, thus it is usually called the 
local preconditioning technique. Many au- 
thors, notably the schools of Merkle [7], Van 
Leer [8]. and Turkel [9]. have made signifi- 
cant contributions in this area. For the sec- 
ond point, the inaccuracies in the upwind 
schemes are primarily due to the incorrect 
scaling of the dissipation term as M — > 0. In 
fact, the dissipation turns out to be scaled 
by the sound speed at low Mach numbers, 
thus yielding excessive numerical dissipation. 
This immediately suggests that numerical 
fluxes need to be modified to correct this sit- 
uation. 

The preconditioning is done to essentially 
alter the eigenvalues of the hyperbolic sys- 
tems so that the wave speeds become more 
or less equalized. To see this, we define a 
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condition number k as the largest ratio of 
eigenvalues, 

j U | 0 

k — — : — : > oc, as ■{/ — > 0, and a held fixed. 

\u\ 

(]5) 

Clearly there is a large disparity of wave 
speeds and as a result this has been identi- 
fied as the source of slow (or no) convergence 
as \v\ — > 0. We bring the attention that the 
limiting taken in Eq. (15) is fundamentally 
different from that by holding u fixed and 
taking a — > oc. The former is more general, 
allowing low speed in a compressible medium 
(where a is finite), and it is especially use- 
ful for dealing with situations with disparate 
speeds of sound in the flow. 

Now if the numerical system is modified in 
such a way that it has a corresponding speed 
of sound, a. which behaves like |u| as it ap- 
proaches to zero, then the condition number 
will remain order of unity. That is. 


1 + a 

t 0( 1 ). while a — > 0, a 

S \ll\ - 

-4 0 

M 




(16) 


Hence, condition number remains of order 
unity at low speeds. The numerical dissipa- 
tion based on this new speed of sound now 
scales with the local speed |t/|. instead of lo- 
cal speed of sound a. To achieve this goal, the 
trick is to manipulate the hyperbolic system 
with the premultiplying matrix. Therefore, 
the system is re-scaled. The numerical speed 
of sound comes in handy for achieving the 
purpose. I Niiig t he time-derivative premult i- 
plying matrix proposed by Weiss and Smith 
[10], the time dependent governing (Euler or 
Navier-Stokes) equations are cast in the fol- 
lowing form: 


mr of da 

1 ~sr + TIT + % - °- 


where W is a vector of primitive variables, 
(P'U.v.T) 1 arid all other variables are stan- 
dard. The preconditioning matrix takes the 


form: 


(") + n T 

uc-) + 7 y; 

r(0 + ~) 
//(0 + 77 r )- 


U p 

pu pc 


P((\ - j) 


e = ?<57J-"' ,l!,) 

Ml = min( 1. max(A/ 2 , A/ 2 ,)). (20) 

The cutoff parameter M co is introduced 
to prevent a singularity at stagnation point. 
It is a user-specified parameter. Unfortu- 
nately it does have some effects on the solu- 
tion in some situations, as will be displayed 
later (in Fig. 4). The effect of M <( , gener- 
ally is minor, but could be of significance in 
some situations. A pressure difference term, 
as suggested by Weiss [11]. can be added to 
enhance robustness near a stagnation point. 
The reference quantity M * is also bounded 
from above to unity as local M exceeds one. 
As a result, the eigenvalues of the flux Jaco- 
bian of F with respect to IU. i.e., V~ ] dF/d\\ 
are //, and 


1 + MI y ( 1 - Mj) 2 M 2 + \M 2 


where M = uja is the unsealed Mach num- 
ber. Several remarks can be made concerning 
the eigenvalues of the preconditioned system. 
Firstlv, we have a bound for the coefficient. 


9 9 — 


Secondly, t he speed of sound is now re-scaled 
by the scaling factor f(M; 4/*). Thus, a new 
speed of sound can be defined. 


a =f{M:M m )a. 


f(M:AL 


- M 2 ) 2 M 2 + 1 M 2 

1 + M 2 


- (24) 
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It is reminded that a in the above equation 
be given in Eq. (12). The scaling factor is 
also bounded. 


1 >f> 


|4/| if 1 >> M 2 » Ml 
2 M co if 1 >> Ml » 


(25) 

That is. the scaling factor / is bounded from 
below by the larger of the local and cut-off 
Mach numbers. Clearly, the variable «, intro- 
duced as above, has been utilized for nothing 
but numerical purposes. Hence it fits in the 
spirit of numerical speed of sound. 

Now equipped with the newly defined nu- 
merical speed of sound, we can readily incor- 
porate it into the formulation of the AUSM- 
family schemes. In this paper, we will il- 
lustrate the concept within the framework of 
A l SM + only: similar procedures can be also 
implemented in the other schemes. Let us 
use the above scaling factor to define the nu- 
merical speed of sound. ii \/ 2 = a{U r U /+1 ), 
at the interface 1/2. The mass flux of the 
ACSM+ scheme in Eq. (10) now can be 
recast by inserting this numerical speed of 
sound. 


ni 1/2— — y- [ M I / 2 ( Pj + Pj + 1 ) + I M 1 / 2 1 ( pj — Pj + 1 ) ] < 

(26) 

where the interface Mach number is given by. 

M l / 2 = jj( Mj ) + ). (27) 

and the arguments in ,Vf ^ have been sub- 
stituted with 

V, = - } [([ + M':).\l J + (\ -A/?).V/,- +1 ], (28) 

A/, + i = l[d t M?)M j+x +(i-Af*)Aij]. (29) 

The tilde Mach numbers. (M r Mj+ 1 ), are 
now defined with respect to the numerical 
speed of sound a, i.e., 

A/ = t. (30) 


This is the scaled Mach number, which will 
revert to the local physical Mach number at 
supersonic speeds. The step taken in Eqs. 
(28)-(29) is necessary to eliminate unphysi- 
cal numerical diffusion present in the pressure 
flux splitting and to enhance the robustness. 
Test cases have shown that the definitions in 
Eqs. (28)-(29) are not strictly necessary for 
the mass flux definition; the simpler alterna- 
tives Mj — Alj and M,+i = Mj+ \ can be 
used as well. These result in a slightly more 
dissipative scheme for low Mach number cal- 
culations. 

For identification purpose, we now call 
this extended method, AUSM + a, to high- 
light the role of t he numerical speed of sound 
a. 

In Fig. 4. we demonstrate that the nu- 
merical dissipation is greatly reduced even if 
the scaled numerical sound speed is included 
in the Van Leer scheme for solving the con- 
tact problem. The error will begin to reach 
the end of the computation domain in 15 
steps with the original scheme, but still re- 
mains inside after 180 steps (not shown in 
the figure) with the scaled numerical sound 
speed. Hence, confirming a much smaller dis- 
sipation as evident in Fig. 4. 
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Fig. I Stationary contact discontinuity, showing 
the effectiveness of using the scaled numerical 
speed of sound in the Van Leer Flux. 

Another interesting result is also revealed 
in Fig. 4. The plot on the right demonstrates 
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the dissipative effect of the cutoff Mach num- 
ber M co on the spreading of the stationary 
contact. The smaller the M co is set. the 
sharper the contact discontinuity gets re- 
sults being closer to the exact solution de- 
noted by solid line, indicating that the nu- 
merical dissipation increases with the value 
of M c0 . 

To further improve the convergence in the 
low speed range, it is found beneficial in Ed- 
wards and Liou [12] that a pressure diffusion 
term rn p be included in the mass flux. 

m 1/2 = Eq. (26)+ mp . (31) 

We can write rn p in the following general 
form. 

- M?)±M( Pj - p, +1 ). (32) 

where 

= [.%. a ,(,V ; ) - M+pj,) 

— -Uj. ;) )f .V/j+i ) + A/j+i )j . (:!-!) 

and the function 'V in the denominator can 
take several forms. Based on the mass flux 
of the AUSMDV scheme, Edwards and Liou 
[12] derived the pressure diffusion term for 
low Mach numbers, 

V = Ml\%- + ^±i]. (34) 

Pj P.,+ 1 J 

In this paper we suggest another formula, 

V = — - — [ + Pj + i )(P:, + P.,+i ) 

PjPj+i 1 

-(Pj - Pj + 1 ){Pj - Pj+l )] • (35) 

This is obtained by also starting with the 
AUSMDV mass flux and by imposing a con- 
dition that its weighting factors be bounded 
(see Appendix A). Furthermore, the last term 
can be omitted to guarantee that P be posi- 
tive. 


P = 


1 


PjPj+i 


( Pj + Pj + 1 ) ( Pj + Pj+l ) , ( 36 ) 


which looks similar to the first expression. 

To see how the numerical speed of sound 
derived for the preconditioned system affects 
the discretization, it is useful to examine the 
mass flux in the limit of an incompressible 
flow, where p is constant and a .\/ 2 (but not 
a i/ 2 ) approaches infinity. Using the defini- 
tions in Eqs. (26)- (32) and taking the limit 
Ml tending to zero in Eqs. (28)-(29), we find 
that 


in 


1 / 2 “ P u \/2 + 

A 


\/"l/2 + 4 ' 


T 2 


V? 


-{pi - Pi+ 1) 


x^(\M,\-\) 2 + 3(Mf- l) 2 ]. (37) 


where 


"1/2 — y + ".,+1 


(38) 


and 


Mi = 


"1/2 


U'U + 4V » J 


39) 


The reference velocity V* is defined analo- 
gously with Ecj. (20): 


U J - max( | U|. U.‘, ). (40) 


This mass flux formula involves only 
the velocity field, pressure, and a constant 
density and is similar to that utilized in 
incompressible-flow discretizations on non- 
staggered grids. Note that the physical sound 
speed completely disappears from the formu- 
lation in the incompressible limit. Advec- 
tive upwind influences are present in the mo- 
mentum equations and the (decoupled) en- 
ergy equation through the switching process 
shown in Eq. (2). but only the pressure- 
diffusion term provides dissipation for the 
continuity equation. 

Moreover, the idea of using a can be read- 
ily applied to solve multiphase flow problems 
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in which both the density and Mach num- 
bers can vary drastically with phase changes. 
More interestingly, the flow can go superson- 
ically even though the fluid speed is low, due 
to a steep decrease in speed of sound. The 
reader is referred to [13] in this conference. 

In what follows we summarize the steps 
involved for adding* the numerical speed of 
sound in t he Al T SM + . These steps, although 
motivated for low- Mach number flows, are 
also valid for high Mach-numher ones. The 
modifications detailed above also work for the 
pressure flux splitting (at least for perfect 
gases), serving to remove 1 unphysical sources 
of numerical dissipation in the momentum 
equations. [12] Further modifications are nec- 
essary for real fluid applications, and the 
reader is again referred to [13] for details. 

1. Fse the numerical speed of sound h\/ - 2 
(23) to define Mj and Mj+\ appearing 
in Eqs. (28)-(29). 

2. Replace the left- and right-state Mach 
number definitions by Mj and A/;+i. 

3. Construct the mass flux as usual by us- 
ing as input the (Mj.Mj+ j) defined in 
the previous step. 

1. If desired, add contribution from pres- 
sure diffusion by using Eq. (32). 

"i. Complete evaluation of other fluxes. 

This is all there is to it , involving a kind of 
pre- and post -processing of the usual AI SM- 
family schemes, steps 1-2 and 4 respectively. 
It is a matter of adding only a few more lines 
to the original A! SM + code. It is believed 
that other low-diffusion hybrid schemes can 
also be extended in a similar manner. We 
shall denote the scheme with t he pressure dif- 
fusion term AUSM+ -ap. 

We now make remarks on the precondi- 
tioning matrix T. We have used the Weiss- 
Smith T to arrive at the scaling function 


f(M:AL) in Eq. (24). Other precondition- 
ers [8,9] can be used as well. The procedure 
for extension will be precisely the same since 
all one needs is the eigenvalues of the pre- 
conditioned hyperbolic system. Thus, a new 
numerical speed of sound a can be expressed 
in terms of the scaling function /(A/ : M x ). 
However, no significant effect on the solution 
is anticipated because all these precondition- 
ers yield more or less the same behavior in the 
limits of Al — > 0 and 1. Unless at low speed 
(say M oo < 0.3), it was found in our calcula- 
tions not necessary to include the precondi- 
tioning matrix in solving the governing equa- 
tions. In other words, the scaling function 
can be incorporated alone, as in Eqs. (26)- 
(29), in the numerical flux and improvements 
in accuracy and convergence can be realized. 

3. Results and Discussion 

In this section, we will present 2D and 3D 
Navier-Stokes solutions for turbulent flows 
over various geometries. The scheme pro- 
posed in this paper was implemented in the 
OVERFLOW code supplied by Bulling et 
al. [14]. The turbulence eddy viscosity in 
all cases presented herein was calculated ac- 
cording to the Spalart-Allmaras one equation 
model [Id]. All of the results presented be- 
low were obtained using an implicit scheme. 
The LHS operator was approximated with 
the standard central difference scheme plus 
appropriate artificial damping terms, (even 
though the RHS residual operator was rep- 
resented with an upwind scheme !). it was 
then further factored and diagonalized in 
each space dimension. 

The flux in the RHS operator was con- 
structed with a third-order accurate inter- 
polation for the primitive variables, together 
with limiter used by Koren [19]. The cut- 
off' Mach number in Eq. (20) is given by 
Ml = Ml/ 4. However, the OVERFLOW 
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code has a parameter controlling the use of 
the preconditioning matrix T. We kept the 
default value to be 3 < 1 under which 
condition T was activated. 

In this paper we will demonstrate the 
effectiveness of using the numerical speed 
of sound in calculating flows at all speeds, 
specifically focusing on two issues: ( 1 ) con- 
vergence rate and (2) accuracy. We will 
show that the convergence rate is improved 
for the entire flow speed regime and the cal- 
culated solutions are in excellent agreement 
with data. 

Shuttle External Tank 

This is an axisymmetric Shuttle external 
tank geometry with a sharp nose and blunt- 
base. downstream of which a significant sepa- 
ration zone is created, see Fig. 8. One of the 


grid lines conforms to the body and grows 
outward and a plane consists of 88 x GO grid 
points. Shown also are the meshes clustered 
around some key regions, one of which is in 
the middle to resolve a tiny notch (not visi- 
ble to the scale). The free stream Reynolds 
number was fixed at 10.000. But a fully tur- 
bulent flow from the nose was assumed in 
order to pose a more stringent condition for 
the assessment of convergence behavior. We 
have tested conditions from low Mach, tran- 
sonic, to supersonic flows. Several schemes 
were considered, consisting of the standard 
AI/SM+ , ATSM+ -a, and AI SM+ -ap. with 
and without the Weiss-Smith preconditioner. 
In all calculations for this problem, we made 
200 steps for each of two coarser grids prior 
to the finest grid, on which 3000 more steps 
were continued unless noted otherwise. 


Table 1: Summary of convergence behavior due to various schemes for the shuttle external 
tank. ft( lX = 10000. 


V, 

Scheme 

5 

£ 

Cu 

0 

X 

P recon d. 


AUSM+ 

v/ 

Diverg. 

0.01 

AUSM+ -a 

No converg. 

y/ 


AUSM+ -a + m p 

No converg. 

yj 


AUSM+ 

yj 

y/ 

0.80 

AUSM+ -a 

V 

y/ 


AUSM+ -a + m p 

c 

y/ 


AISM + 

yj 

yj 

2.00 

AUSM+ -a 

yj 

sj 


AUSM+ -a + ni r 

y/ 

si 


Table I summarizes the convergence be- 
havior of the above combinations. We ob- 
serve the following: (1) For low Mach num- 
bers (approximately yV/ x < 0.3). it was found 
necessary to use the time-derivative precon- 
ditioner T so that the numerical dissipations 


in both the implicit and explicit operators 
are compatibly scaled. Otherwise, the cal- 
culation either diverged or stagnated. (2) 
For flows at transonic speeds or higher, the 
time-derivative preconditioner, as given in 
Eq. (18). serves no benefits whatsoever, even 
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though the fluid speed is low in the vis- 
cous and the base flow regions. This is un- 
derstandable because the preconditioner ef- 
fects only t he inviscid waves and the informa- 
tion in t he viscous-dominated regions is only 
transmitted via diffusion processes which are 
ably handled by the implicit operator. Vis- 
cous arid grid aspect effects can be included 
in the construction of preconditioner, see for 
example [Hi]. 



Fig. a Convergence history for the shuttle exter- 
nal tank problem. 


In Fig. 5. we display the convergence 
history for various Mach numbers using 
AISM + . without using the preconditioner 
F. The residuals for the low Mach-number 
cases stall after a drop of four orders of mag- 
nitude. These drops in many calculations, al- 
though not especially admirable, would have 
been acceptable. However, a close examina- 
tion of the solution reveals that it is com- 
pletely unacceptable, as shown in Fig. G. It 
appears that then' is a false boundary (ex- 
actly aligned with a grid line) at which infor- 
mation is unable to pass. This phenomenon 
is quite typical in the low Mach-nurnber cal- 
culations using an unmodified compressible 
code, also seen in [12]. Hence, a measure 
of caution should be taken when reading the 
residual history for the low Mach-number 


solutions. 


\ \ 




Fig. 6 Pressure contours for the shuttle external 
tank problem obtained at N=6400 time steps for 
M iyi , = 0.01, using the standard AUSM + . The 
Bottom picture shows a magnified view near the 
nose. 


i c° 
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Fig. 7 Convergence history for the shuttle ex- 
ternal lank problem obtained by the AUSM+ -a 
scheme. 

On the other hand, the convergence his- 
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lories with the use of the numerical speed of 
sound display improvement, over those with- 
out it. as shown in Fig. 7. As noted ear- 
lier in Table I. it is necessary to invoke F for 
the \/ x = 0.0 land 0 .1 cases. The conver- 
gence rates for these two calculations nearly 
coincide with each other, indicating Ma.cli- 
number independence. 


- ■. , .. • • i i i <’*£*// 



Fig. 8 Pressure contours for the shuttle external 
tank problem obtained at N = 1000 time steps for 
il/, x . — 0.01. using AUSM+ -a. The middle pic- 
ture shows a magnified view near the nose and 
the bottom one depicts the separated zone be- 
hind the tank base. 

In Fig. 8, we show the solution at N = 1000 
* However, t he pressure distribution along the wall is 


steps at which the residual has been dropped 
to the level approximately equal to that 
shown in Fig. 6 (N=6400). It is of inter- 
est noticing that the solution now is well be- 
haved and is every bit as good as the solu- 
tion at N=3200. Also the blow-up view near 
the surface reveals smooth profiles of pressure 
contours, unlike the standard AUSM+ which 
has been known to yield unwanted pressure 
oscillations in viscous layers along the trans- 
verse grid lines when the mesh aspect ratio 
is large and flow is essentially parallel to a 
grid line A The separation region in the base 
of the tank is depicted by particle traces, in- 
dicating its size extending about one radius 
downstream. 



2CC0 4000 6000 


M 

Fig. 9 Convergence history for the shuttle exter- 
nal tank problem obtained by the A l SM + ap 
method. 

Finally ihe effect of including the pres- 
sure diffusion term on the solution was inves- 
tigated and the results are given in Fig. 9. 
Again, the preconditioner F must be used for 
the low Maeh-number cases and their conver- 
gence histories are essentially identical, be- 
coming independent of Mach number as the 
Mach number lowers. The pressure contours 
are indistinguishable from those shown in 
Fig. 8 and are t hus not included. 

smooth. 
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Fig. 10 Pressure contours at = 

0.01.0.1,0.8.1.25, and 2.0 (from top to bottom) 
for the shuttle external tank problem, showing 
accordingly different flow features. 

The solutions obtained by using the 


AUSM + a for several A/, x values are given 
in Fig. 10 for comparison. They all ap- 
pear physically correct and are numerically 
well behaved. No discenible differences were 
noticed between those by AUSM + -a and 
AUSM+ -ap. 

Comparing Figs. 5, 7 and 9, we see that 
the convergence rate is improved in the tran- 
sonic ranges by simply using the numerical 
speed of sound alone. For low Mach num- 
ber cases. M x = 0.01,0.1, another order of 
reduction can be obtained by including the 
pressure diffusion term. Also, the use of nu- 
merical speed of sound yields the convergence 
histories that are relatively insensitive to the 
flow speeds. Thus, the validity of the present 
method is confirmed and the goal of having 
a convergence rate that is more or less inde- 
pendent of the free- stream speed is achieved. 

We now summarize major findings from 
the study of this problem: (1) The numer- 
ical speed of sound concept is an effective 
means of extending AUSM-tvpe discretiza- 
tions to solve low Mach number flows in an 
accurate and efficient manner. (2) Since the 
numerical speed of sound is reduced with the 
flow speed, the numerical dissipation changes 
accordingly, and a compatible implicit op- 
erator (one that includes the precondition- 
ing matrix) must be used. (3) For moderate 
Mach numbers and beyond, it is not neces- 
sary to use T. (4) Incorporation of the nu- 
merical speed of sound, as described in steps 
(26)- (29). helps remove pressure oscillations 
in the viscous layers. 

ONERA M6 Wing 

The next problem is the ONERA M6 wing 
with the free stream conditions A/ x = 0.84, 
and /?c x = 18.2 x 10 (> . under various angles 
of attack. The computation domain consists 
of 269 x 35 x 67 grid points. For this case, 
the preconditioning matrix V was automati- 
cally turned off in the code since the control- 
ling factor 3A/,*, exceeds unity. However, the 
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numerical speed of sound a was active with 
Al co = My:./ 2. The pressure contours on the 
wing surfaces are displayed in Fig. 11, show- 
ing the well-known A-sliock pattern, thus ap- 
pearance of two shocks roughly in two thirds 
of the wing span. The detailed comparison of 
surface pressure distributions are shown for 
two spanwise sections. Figs. 12 and 13 re- 
spectively for 44% and 65%, for four angles 
of attack. The computed results are in very 
good agreement with the data [17], especially 
in capturing of the shock locations. 

The convergence histories are presented 
in Fig. 11 for two different angles of attack. 
They show a continuing decrease by about 
five orders of magnitude from the largest val- 
ues. at nearly the same rates. 

Wingbody 

Turbulent flows over a wingbody config- 
uration are calculated and their results are 
now discussed. The geometry is shown in 
Fig. 15. where the sting is included in the 
calculation. The computation domain is grid- 
ded using the chimera overset grid technique 
and entire grid is composed of seven grids. 
The flow conditions are: M iyi — 0.8, o = 2°, 
and /fe lX , = 0.167 x 10T Figure 16 depicts 
the pressure coefficients at various spanwise 
locations. The computed results are in excel- 
lent agreement with the measured data [18]. 
Moreover, the pressure coefficients along the 
body, shown in Fig. 17. exhibits a similar 
level of excellent agreement with the data, 
even in the wing root region where a sharp 
variation is encountered. 

Finally, Fig. 18 displays a well-behaved 
convergence history, reducing the residual er- 
ror by more than five (5) orders of magnitude 
in 800 steps. 

4. Concluding Remarks 

In this paper we have introduced the no- 
tion of "numerical speed of sound'", which 


turns out to be a very useful variable that can 
be employed to construct an upwind scheme 
to satisfy certain properties, most notably 
the exact capturing of contact and shock dis- 
continuities for ID problems. This numer- 
ical speed of sound is explicitly utilized in 
the construction of AUSM + and AUSMDV 
schemes. However, the idea can be inserted 
in the other upwind schemes. One example 
is in the Roe flux splitting where an aver- 
aged speed of sound, among several other av- 
eraged variables, is automatically required for 
the process. We have shown that other cele- 
brated schemes, which in standard form ex- 
hibit intermediate points, can now be made 
to capture a shock exactly. 

More important ly, the concept of numer- 
ical speed of sound, being meaningful only in 
the numerical sense, can be extended to effec- 
tively deal with flows at low speeds. The crux 
is that a scaling factor varying with speed (or 
Mach number) is introduced. As a result, the 
numerical dissipation is decreased with the 
flow speed. Hence the convergence rate is en- 
hanced, not only at low speeds, but also at 
high speeds as well. Additionally, the solu- 
tion accuracy is improved. 

The effectiveness of implementing the nu- 
merical speed of sound in the AUSM + scheme 
has been demonstrated. Solutions of com- 
plex turbulent flows were obtained for com- 
plicated geometries, meshed with the over- 
set grid technique, using the OVERFLOW 
code. We have presented convergence his- 
tories demonstrating significant improvement 
over the previous scheme. The pressure dis- 
tributions are' in excellent agreement with 
available data for the ONER A M6 wing and a 
wingbody configuration, further proving the 
reliability of the new scheme, A1 SM + a. 
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Appendix A 

In this section we give a formal derivation 
of the function T? that appears in the pressure 
diffusion term in mass flux. We begin with 
the AUSMDV flux written in terms of the 
numerical speed of sound. 

m i/ 2 — h \ /2 ) ( M j ) + ) 

+ Pj+ l*' + l ))], (1) 

where 
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These can be rewritten as 


PM + = tf:(I 2p-±p), 


Pj + = ^q(12p+ -Vt 


together with 


As <1 result , we should have defined u, ,:t in the 


following way. 


2 D^ P 


M 2 h 


P.i + i 




Y,P = Pi + PJ + 1 ■ A /> = 


/h+i - Pj 


Noticing that D has the dimension of p/p. 
Substituting these formulas into Eq. (1), we 
denote the coefficient for Ap term by 


I<p = 


-AM. 


where we have A.V( given in Eq. (33). As 
M r M ;+ 1 — ► 0. we get 

A.Ad = (- + T). (7) 

Hence, if D — 0(1). then 

K p = 0(t/ 1/2 ), as u i /2 < 1. (8) 

That is. the pressure diffusion term dimin- 
ishes with the flow speed arid it clearly is in- 
sufficient to provide adequate contribution of 
pressure to the mass flux at low speeds. This 
can be remedied by properly rescaling I\ p as 


= wrrM- 


What is left is to determine the expression of 
D. First, we remark that the condition for 
capturing a stationary contact discontinuity 


requires 


Pj*' ~ Pj+i l 


Monotonicity constraints on the and 
split Mach numbers in Eq. (1) implies condi- 
tions on a.’* (see [4] for example), 


“ - 1 + lb' v 7 

A more stringent condition fulfilling the 
above inequalities may be 


0 < + , 


This gives the relation 


s TT7v “ j > 


V = M*D = 


PjP.i+i 


•V ?(]>.! + P.i+ 1 ) ( Pj + P.i+i ) 


-(P.i - Pj+i )(pj ~ Pj+ 1 )] •( f r > ) 


This completes the derivation of D. 
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Fig. 1 1 Pressure contours on the ONERA Mft wing at a — 3°\ A/ x = 0.84, and iZc, x . = 18.2 x 10 H , 
showing the A-shock pattern near the wing tip. 


_ ' 5 , 




-1.5 

' -C 



_ - Exst 

CFO 

(\= 0.04 de$ 

-1.0 

- 0.5 : 

: ^ 



- 0.5 

c.c ■ 

i 


'*' ; e . ; 

CL 

CJ 

0.0 

C .5 ■ 




0.5 





1 .0 

. u 

0.0 

0.4 

x/c 

0.8 

" t; 




1 s 




"" Exot 

f . S_J 

.0 



a - 4.06 detj} 

- 1 .0 

- 0.5 

C.C 

- 

fr . -♦ " ■■■■' 

. ! 
: 'K.. ! 

- 0.5 

CL 

0.0 

C .5 ■ 




0.5 

- r\ , 




1.0 

. C/ 

o 

o 

0.4 

C.S 




0.0 


).0 


0.4 


x/c 


r Ex 3* i 

— c r o 

a— 3.06 de<j 


0.8 


C Exd* j 

- c r o 

n- 5.06 de§ 


1 C-.N 


0.4 


0.8 


x/c 


x/c 


Fig. 12 Pressure distribution at the span wise section. 44%. of the ONERA MG wing for various angles 
of attack. ,\7 X = 0.X4. 7fc x = 1X.2 X 10'\ 
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x/c x/c 

Fig. 13 Pressure distribution at the spat) wise section. 65 ( X. of the ON ERA M6 wing for various angles 
of attack. M x = 0.84. R< x = IX. 2 x 10 6 . 



Fig. 13 Convergence history for the ONERA MG wing problem at two angles of attack. 4/ x = 
0.84. fle*. = 18.2 x 10'\ 
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Fig. 15 Geometry of the wingbody problem. 
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Fig. 16 (rout'd) 



x/c x/c 

Fig. 17 Pressure distribution at angles along the body. A/ x , = 0.8. n = *2'\ R r lX , - 0.167 X 10 <s . 
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Fig. 17 (coat'd) 
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Fig. IK Convergence history for the wingbody problem. 
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